Introduction

Macroscopic mechanical behaviour of saturated sand is significantly influenced by several internal and external factors. The internal factors include particle strength [33, 68, 83], particle size [16, 85, 87], particle size distribution [17, 59, 75] and particle shape [20, 37, 68, 79, 80], while the external factors include confining pressure [7, 59, 76], stress path [66, 78], loading mode [10] and drainage conditions [11]. In nature, sands are not clean since a certain amount of fines are added with particle size finer than 74 μm. The presence of fines is considered to change the internal parameters. It can greatly affect the strength, dilatancy and critical state behaviour of these sand–fines mixtures in both drained and undrained conditions [47, 52, 54, 56, 57, 63]. Since many interfacial factors control the mechanical behaviour of silty sand, the different or even contradicting results reported in the literature indicate that the influence of fines remains uncertain [53]. As a result, very diverse views exist as to whether the effect of fines on the shear strength of silt–sand mixtures is negative or positive [80]. Murthy et al. [45] performed a series of undrained monotonic triaxial compression tests on Ottawa sand mixed with crush silica. They concluded that an increase in non-plastic fines content of sand leads to an increase in the critical state friction angle. Similar conclusions were reported by Wei and Yang [69], Chang et al. [6], Pitman et al. [49] and Kuerbris et al. [31]. On the other hand, based on their triaxial tests on quartz sand mixed with kaolin, Georgiannou et al. [21] concluded that the degree of collapsibility of mixed soil increased with the addition of a small amount of fines.

The degradation in shear strength of soil because of fines can be explained by the binary packing representing the structure of coarse and fine mixtures proposed by Lade and Yamamuro [36]. Nevertheless, their round-to-round model cannot explain the gain in shear strength observed in other studies when fines are added. Recently, Yang and Wei [80] used an angular-to-angular model to interpret the increase in the stability of the soil structure when increasing fines content. They compared undrained shear behaviour of different mixtures of host sand (Toyoura/Fujian sands) and fines (crushed silica/glass beads) and quantitatively studied the influence of particle shape. They showed that the undrained shear behaviour and collapsibility of a sand–silt mixture are closely related to the shape of its constituent particles: a mixed soil containing rounded fines tends to exhibit a higher susceptibility to liquefaction collapse than a mixed soil containing the same percentage of angular fines. Several other experimental and numerical studies also support the crucial influence of particle shape on the behaviour of granular materials [30, 44, 68, 73, 79, 86, 88], including a recent 2D DEM study by Gong et al. [23], who explored the clumping particle technique to investigate the effect of particle shape on the shear behaviour of silty sands.

Undoubtedly, the particle shape has a significant impact on the complex mechanical behaviour of sand–silt mixture, which is essentially due to the heterogeneous and discrete nature of the soil that leads to the complicated internal structure under shearing. However, there are multiple definitions of shape parameters found in the literature such as sphericity, aspect ratio, convexity, overall regularity, elongation and slenderness [10, 11, 77]. The shape evaluation based on these parameters is not only complicated but also challenging to achieve in grain-scale modelling to evaluate the internal texture of the soil. For simplicity, rolling resistance has been incorporated into DEM formulations to account for a rotational frictional torque caused by particle shapes. The rolling resistance represents not only the influence of particle shape but also micro-slip and friction on the contact surface, viscous hysteresis and surface adhesion [1]. Therefore, as realised by many researchers [1, 20, 27, 43], rolling resistance plays an important role, particularly in DEM simulations of circular or spherical particles. However, the effect of rolling resistance on the shearing behaviour of silty materials has not been thoroughly analysed yet.

This paper examines the effect of particle rolling resistance on the shearing behaviour of silty sand under both drained and undrained conditions using 3D DEM simulations. The rolling resistance of fines and coarse particles along with fines contents is varied in triaxial shear simulations to enable a better understanding of the mechanical responses of silty sand at different states during triaxial loading.

Numerical setup

A series of triaxial tests were numerically conducted using DEM to study the effect of particle rolling resistance on the shearing behaviour of silty sand mixtures. DEM has been widely used to investigate the mechanical behaviour of granular materials since its first introduction [15]. It is a precious tool to capture local responses and micromechanical details, all of which significantly influence the macro-behaviour of the granular materials. Several DEM models have been proposed in the literature for different purposes [14]. In this study, the Hertz–Mindlin contact model was adopted to describe the force-dependent contact stiffness [24, 25]. The elastic–plastic spring–dashpot (EPSD) model was employed to provide a stable resistance torque in quasi-static condition [1]. The numerical details of the DEM model can be found in Appendix A.

The specimen used in DEM triaxial simulations was represented in a cubic RVE (approximately \(3\;{\text{mm}} \times 3\;{\text{mm}} \times 3\;{\text{mm}}\)) by three pairs of stress or strain-controlled rigid and frictionless boundaries. The specimen size is kept small to reduce the computational time and is large enough to capture the average material behaviour while minimising the boundary effect as it is larger than 18 \(\times\) mean particle size \(D_{50}\) and has the number of particles varying between 10,000 and 70,000 [12, 25, 32, 50]. The basic parameters for DEM simulations are summarised in Table 1. The selection of rolling resistance coefficients is explained in Appendix B.

Table 1 Input parameter for DEM model

Figure 1 shows the particle size distribution (PSD curve) used in this study. The PSD of host sand was selected to fit the central portion of the PSD of Nevada sand. The largest particles (3.5%) were omitted because they have less contribution to the overall behaviour of material [32]. The PSD of fines particles is in the range of coarse silt. Therefore, the sand–fines mixture has a well-graded PSD in overall and the size ratio between the mean particle size \(D_{50}\) of fines and coarse particles is around 3.5, and the ratio between the largest and the smallest elements is around 6. Parallel computing is employed to reduce computational time due to the disparity of the particle sizes when adding fines.

Fig. 1
figure 1

a Particle distribution for sand and fines; b DEM assemblies for different fines contents

The specimens were prepared by random distribution and rearrangement of DEM particles in the RVE. First, a number of spheres were randomly generated in 3D cubical space which is 5% bigger than the expected size of the sample with the void ratio of 0.4 without checking the overlap among particles. The overlap will gradually smooth out when particles rearrange. In this period, the motion of particles is limited within 1% of their radius and the friction coefficient is set to infinity to overcome explosion due to high overlap at the beginning. The initial state is reached by isotropically confining the specimen to 10 kPa. During this process, different particle frictions were altered to generate the samples with different densities [25, 71]. In general, a low value of friction coefficient leads to a dense specimen and vice versa.

Particles shape plays an important role in the behaviour of granular materials; however, it is also one of the most challenging properties to model accurately. Therefore, in this study, rolling friction is employed to introduce ‘shape-like’ behaviour to DEM simulations simply. A series of monotonic drained and undrained triaxial tests were conducted for different fines contents, initial void ratios and confining pressures to estimate the effect of rolling resistance as well as shape factor, correlated with fines content, on the behaviour of silty sand. These tests are divided into four different cases (A to D), and the material properties corresponding to each case are shown in Table 1. In each simulation case, the percentage of the fine (\(f\)) was varied from 0 to 50% to observe the effect of fines content in combination with rolling resistance on the behaviour of sand–silt mixture. The drained tests for each case were simulated with five different confining pressures \(( p_{0} )\) of 100, 200, 300, 400 and 500 kPa, and different initial void ratios (\(e_{0}\)) ranging from loose to dense states. The confining pressures are relatively low to avoid underestimating the interlocking effect.

Effect of particle rolling resistance on drained behaviour of silty sand

Depending on whether the initial state is dense or loose, the shear behaviour of soil in the drained triaxial tests can be divided into two types, dilatant and contractant, some of which are shown in Fig. 2. The loose sample exhibits contractant behaviour where the volumetric strain continuously decreases, and the shear stress continuously increases during the shearing process until reaching its critical state where the sample continues to shear at constant effective stress and volumetric strain. The dense sample shows similar behaviour at the beginning of shearing with slight volumetric contraction. After this, volume dilation takes place until the peak state is reached. The shear stress then softens towards the critical state. The drained behaviour of soil is, therefore, characterised by the critical state, peak state and zero-dilatancy state (also known as the characteristic state).

Fig. 2
figure 2

Shear and volumetric behaviour for different fines contents for medium dense (continuous lines) and medium loose samples (dash lines) at confining pressure of 100 kPa

The results shown in Fig. 2 highlight that the particle rolling resistance affects the drained shearing behaviour of silty sand strongly. As shown in Case D, samples with higher fines contents show higher dilatancy, higher peak strength and higher residual strength if fines particles have higher rolling resistance than coarse particles. Otherwise, if fines particles have lower rolling resistance than coarse particles, increasing fines content reduces the dilatancy, peak and residual strength of soil as can be seen in Case B. Moreover, data for Case A and Case C show that when coarse and fine particles have the same rolling resistance, adding fines causes less effect on the overall behaviour of well-graded silty sands. In other words, the trends in both shear resistance and volumetric behaviour in such cases (A and C) remain unchanged, despite noticeable effects on the shear strength and volumetric strain. Undoubtedly depending on the correlation between rolling resistance of coarse and fines particles and also the percentage of fines in the sand (Cases B and D), the critical state, peak state and zero-dilatancy state along with the trends in both shear resistance and volumetric behaviour vary significantly. The detailed effect of rolling resistance which represents the effect of particle shape on the drained behaviour of silty sands will be discussed in the following sections.

Critical state behaviour

The critical state has been defined as the state at which the soil continues to shear at constant stress and constant void ratio [55]. The locus of critical state points on the deviatoric stress (\(q\)) and the mean effective stress (\(p^{\prime}\)) plane or the mean effective stress and void ratio (\(e\)) plane forms the critical state line. The critical state line is generally considered to be unique and independent of stress path, sample preparation method and initial density. Therefore, the critical state line provides a reference state from which the state parameter \(\psi\) and other important parameters of soil are derived. The critical state line in the \(p^{\prime} - q\) shear plane is specified by a straight line passing through the stress origin with a slope \(M\):

$$p^{\prime} = Mq$$
(1)

While in the volumetric plane or \(e - p^{\prime}\) plane, the critical state line is represented by the power-law function [39]:

$$e = e_{\varGamma } - \lambda_{c} \left( {\frac{{p^{\prime}}}{{P_{A} }}} \right)^{\xi }$$
(2)

where \(e_{\varGamma }\) is the maximum void ratio, \(\lambda_{c}\) is the slope of critical state in volumetric plane, and \(\xi\) is the pressure exponent. In this study, \(\xi\) is taken the value of 0.7, which is the best-fitted value to the test data.

The values of \(e_{\varGamma }\) and \(\lambda_{c}\) are unsurprisingly found to be strongly affected by fines content, as shown in Figs. 3 and 4, respectively. Increasing fines content up to a threshold value (also known as transition fines content—TFC) will reduce the maximum void ratio \(e_{\varGamma }\). However,\(e_{\varGamma }\) will increase with the fines top up if the amount of fines is higher than the transition fines content. These changes in the compressibility of the silty sand mixture are consistent with reports published in the literature [47]. This can be explained with binary packing consideration in which fine particles are supposed to fill up the void space among coarse particles. In fact, the binary packing mixtures exhibit a shaper transition where the TFC value is ranging from 20 to 50% depending on the particle diameter ratio as reported by Rahman [53]. When the soil is well graded, the transition fines content is within 20–30% and the curvature is relatively low because the distinction between the fine and coarse fractions is comparatively small and the void space by coarse particles is smoothly filled up with fines [42, 45, 48, 81]. In this study, as we used well-graded PSD, the TFC can similarly be observed at around 30% of fines. Interestingly, TFC is less likely to depend on particle rolling resistance. This can be primarily attributed to the fact that the fulfilment of fines into the void space is not primarily controlled by particle shape but by particle volume or particle size distribution.

Fig. 3
figure 3

Critical state lines in the volumetric plane with different fines contents and rolling resistances

Fig. 4
figure 4

Comparison of critical state parameters in volumetric space: a maximum void ratio; b the slope of the critical state line on the volumetric plane

Although causing less influence on TFC, the rolling resistance strongly affects the compressibility of sand–silt mixture where it widely alters the maximum critical void ratio \(e_{\varGamma }\) and the slope of the volumetric critical state line \(\lambda_{c}\). This can be explained by the change in compressibility of soil with the presence of fines contents causing changes in the interaction network. When fines are added into the mixture, the number of mutual contacts among coarse and fine particle including coarse–coarse, coarse–fine and fine–fine contacts significantly changes. Regardless of particle shape, increasing fines percentage coincides with the growth of coarse–fines and fines–fines contacts and the reduction of coarse–coarse contact as shown in Fig. 5. Similar behaviour can be observed with the clump particles in the study of Gong et al. [23]. Besides, rolling resistance has a minor effect on the proportion of contact types in the mixture, which is clearly shown in Fig. 5a. Also, in this figure, the TFC can be found close to the intersection of coarse–coarse and fines–fines lines. Nevertheless, when the strength of these contacts vigorously changes along with particle rolling resistance, it will lead to the change in the overall strength of the mixture correspondingly. Results for Case A and Case C in Fig. 4 show that as rolling resistance of mixture simultaneously increases, the \(e_{\varGamma }\) − \(f\) line shifts up in parallel from black line to red line. This indicates that rolling resistance enhances the contact strength among particles, and therefore, maximum critical void ratio \(e_{\varGamma }\) is higher regardless of the fines content. However, if the rolling resistance of coarse particles is kept constant and rolling resistance of fines particles is reduced, \(e_{\varGamma }\) will considerably decrease as fines increases. This degradation is associated with the growth of coarse–fines and fines–fines contacts, which are weakened as rolling resistance drops. Inversely, if rolling resistance of fines particle is relatively increased compared to coarse particle, \(e_{\varGamma }\) will increase as fines increases. However, the variation of \(e_{\varGamma }\) associated with the high rolling resistance of fine particles is less significant than that with low rolling resistance fines particles. This is due to the complex correlation between rolling and sliding friction in addition to the variation of volumetric critical state slope \(\lambda_{c}\) as can be seen in Fig. 4b. Rolling may play a dominant mechanism at low rolling resistance; however, at the higher rolling resistance, sliding mechanism mainly controls the shearing process especially at low confining pressure [18]. Wensrich and Katterfeld [70] have also reported similar observation for both spherical particles with rolling resistance and clumped particles.

Fig. 5
figure 5

Coordination numbers at critical state for 100 kPa confining pressure: a Percentage of contact types is independent of rolling resistance; b coordination number for fines–fines contacts; c coordination number for coarse–coarse contacts; and d coordination number for coarse–fines contacts

Based on the available experimental results reported in the literature, the influence of fines on the volumetric critical state slope \(\lambda_{c}\) of a sand–silt mixture can be grouped into two broad views: (1) slope significantly rotates around a pivot point with increasing fines content and (2) slope is more or less parallel at certain stress level [81]. The numerical results reported in this paper inclined towards the second view. Accordingly, the volumetric critical state line for cohesionless particles on well-controlled gradation will slightly rotate with the change in the amount of fines added, and no pivot point has been identified. The critical state line somewhat rotates counter-clockwise when increasing fines if fines percentage is smaller than TFC and moves inversely if fines percentage is higher than TFC. This can be explained by the variation in the compressibility of soil with the mutating in the mixture volume by filling fines particle incorporated with the change in the contact network. When fine particles contribute a larger proportion in the mixture, they play an important part in soil skeleton rather than filling void space; therefore, higher confining pressure is required to reach the critical state. The balance between the contribution to soil skeleton and the occupation of void space by fine particles is dedicated leading to the trivial change in the critical state slope. From Figs. 3 and 4, it can be concluded that fines rolling resistance does not play an important role in rotating the critical state in the volumetric space. High rolling resistance fines more strongly alter \(\lambda_{c}\) rather than \(e_{\varGamma }\), and no pivot point has been found as reported by Been and Jefferies [2] and Bouckovalas et al. [5]. Similar to the interlocking mechanism, increasing the confining pressure will enhance the effect of rolling resistance on the particle rotation. Thus, the contribution of the sliding to the shearing process becomes more and more significant. However, the coupling between sliding and rolling resistances is quite complex as evidenced by the lower \(\lambda_{c}\) in case D than in case C when increasing fines content with higher rolling resistance. This behaviour can be further explored to clear the gap between the volumetric behaviour of clumped and spherical particles displayed in Fig. 22 in Appendix B.

The effects of rolling resistance on the critical state in the stress space are shown in Fig. 6a, b. Figure 6a shows the critical state lines interpolated from the critical stresses in the series of drained tests in Case B. In contrast, Fig. 6b shows the plot of critical state slope (\(M\)) in the shear plane, also known as critical stress ratio, versus fines percentages. It can be seen from Fig. 6a that the shear resistance of sand–fines mixture in Case B notably reduces when adding fines which have lower rolling resistance than that of the original sand. The critical state in Case B tilts in the clockwise direction similar to the observation of Marto et al. [42], Phan et al. [64], Thevanayagam et al. [63]. The shear resistance of the mixture is mainly governed by the particle contact types and the fraction of those contacts within. As shown in Fig. 5, the proportion of different contact types is almost constant and independent of rolling resistance. However, the strength of coarse-to-fines and fines-to-fines contacts significantly increases with the increase in the anti-rotation effect. Thus, it can be seen that higher rolling resistance of fines enhances the shear resistance of the mixture when comparing the results among Cases B, C, and D in Fig. 6b. In reverse to Case B, the critical state slope \(M\) in Case D increases with fines increase. This result is consistent with the literature [6, 45, 69], which reported that increase in fines content (even small amount) enhanced the shear resistance of sand–silt mixture. The substantially different contact properties between coarse and fines strengthen the effect of fines on the mixture shear strength. Fines with low rolling resistance in Case B form an unstable or metastable structure which makes the large particle held slightly apart near contact points and slide to each other. In contrast, fines with high rolling resistance in Case D strengthen the coarse-to-fines and fines-to-fines contacts and, therefore, reduce the instability in the load-bearing structure of sand–silt mixture.

Fig. 6
figure 6

Drained critical state in stress space: a critical state in Case B; b comparison of critical state slopes among cases

Undoubtedly, different rolling properties between coarse and fines particles result in a significant change in shear resistance. However, if fines and coarse particles have a similar anti-rotation level, adding fines has a trivial impact on the shear strength of the mixture. This can be explained by the similarity in contact properties among coarse and fines particles despite the difference in the contribution of contact types in the force chain network shown in Fig. 5. As shown in Fig. 6b, when rolling resistance of fines is equal to that of coarse particles (Case C), adding fines does not noteworthy alter \(M\). \(M\) slightly goes up when increasing fines percentage up to 20% and then slowly decreases for a further increment in fines.

Similarly, when no rolling resistance is assigned to all particles (Case A), \(M\) somewhat decreases until the percentage of fines reaches 30% and then plateaus down. Although the trends of \(M - f\) lines in these two cases are different, the variation is trivial, and these two lines are almost parallel. This is because of the analogous percentage of contact types at the identical fines content, as shown in Fig. 5a. Therefore, it can be concluded that when increasing the rolling resistance of the whole mixture, \(M - f\) line is parallelly shifted up. The current results show that the shear strength in cohesionless soil is strongly governed by particle shape, which is modelled by the plastic rolling moment, in addition to the inter-particle friction angle as previously reported by [19, 43].

In summary, the rolling resistance of fines strongly affects the critical shear and volumetric behaviour of sand–silt mixture, which mimics the influence of particle shape on the mechanical behaviour of sand–silt mixture commonly reported in the literature. They cause changes in particle contact strength along with the variation of the number of contacts by adding fines. Rolling resistance also strongly affects the coordination number of coarse–coarse, coarse–fine and fines–fines contacts. However, the percentage of these contacts in the mixture is not affected by rolling resistance. As a result, simultaneously increasing the rolling resistance of fines and coarse particles will reduce the compressibility of the mixture. When fines and coarse particles have similar rolling resistance, the critical state line in the shear plane does not significantly alter with varying fines content. Nevertheless, a relative difference in rolling resistance between coarse and fines particles strongly affects the CSL in both shear and volumetric planes. In particular, for the same particle size distribution, higher ratios between rolling resistances of fines and coarse contents lead to improvements in shear resistance. However, it is noticed that rolling resistance has a minor effect on altering the value of TFC, which is close to the where coarse–coarse and fines–fines contact percentage lines meet. Also, it does not play an important role in the rotation of the volumetric CSL. The strong rotation of CSL in volumetric plane observed in the literature may not result from particle shape. Furthermore, the notable rotation of critical state in the shear plane shows that inhomogeneous rolling resistance incorporation with sliding friction can be used for a better presentation of particle shape with spherical particle.

Peak state

A peak state is defined as the state where the deviatoric stress is maximum. The axial strain at the peak state corresponds to the maximum rate of dilation, which is referred to as the state where the ratio of the volumetric-strain increment to the axial-strain increment is maximum [4]. The peak phenomenon is effectively suppressed, while the sample volume is contracted rather than dilated. Because fines with different rolling resistances significantly alter the volumetric as well as shear behaviour of sand–silt mixture, dilatancy and hardening responses of the mixtures are strongly affected by rolling resistance of particles.

Results from a series of drained triaxial tests for medium dense samples shown in Fig. 8 exhibit a strong dependence of peak state on the rolling resistance of fine particles, as compared to the confining pressures and fines contents. The peak strength of sand is often represented by the peak friction angle \(\phi_{p}\), which is defined by major and minor principal effective stresses at the peak state \(\sigma_{1p}^{{\prime }} , \sigma_{3p}^{{\prime }}\):

$$\sin \phi_{p} = \frac{{\sigma_{1p}^{\prime } - \sigma_{3p}^{\prime } }}{{\sigma_{1p}^{\prime } + \sigma_{3p}^{\prime } }}$$
(3)

where \(\phi_{p}\) refers to the ‘secant’ angle of shearing obtained by dropping a tangent from the origin on to a single Mohr circle of effective stress [4]. It is evidenced from Fig. 8 that if rolling resistance of fines is relatively smaller than that of host sand, increase in fines percentage will reduce the peak strength of the mixture. This can be attributed to the growing number of round-to-round contacts with low shear strength. In contrast, if rolling resistance of fines is comparatively higher than that of host sand, adding fines is equivalent to increasing the number of strong contacts, and hence, the shear strength. This observation agrees well with the experimental results of Xiao et al. [76, 77] when considering the effect of particle shape on peak state by adding glass bead fines with a different roundness.

The single plot of peak friction angle for different rolling resistances of fines can be found in Fig. 8b–d. There is a notable trend when fines content is higher than TFC. When fines rolling resistance is fairly higher or lower than coarse rolling resistance, the peak friction angle will plateau down for a further increment of fines percentage after TFC, i.e. 6 (b) and (d). Moreover, TFC again plays the role as the transition milestone of peak friction angle if rolling resistance of fines is similar to that of host sand, as shown in Fig. 7c. This can be explained by the strong correlation between fines added and the change in soil structure. The saturation of void space filled up with fines particle will lead to the major occupation of fine–coarse or fine–fine contacts which are controlled by rolling resistance of fines. While it is well known that particle rotation plays a dominant role in the development of shear [56], the peak state after TFC will approach the state of the ultimate fines content (100% fines).

Fig. 7
figure 7

Peak friction angle with different fines rolling resistances: a comparison among 3 different rolling resistances of fines for Cases B, C, and D; b Case D; c Case C; and d Case B

Soil often exhibits hardening and softening responses before and after the peak, respectively. Those behaviours are considered to be strongly dependent on two distinct state variables: state parameter \(\psi\) and stress parameter \(\eta /M\). The state parameter \(\psi\) is the void ratio difference between the current state and critical state at the same mean effective pressure (\(p^{\prime}\)): \(\psi = e - e_{c} \left( {p^{\prime}} \right)\), while the stress ratio \(\eta\) is the ratio between shear stress \(q\) and mean effective stress \(p^{\prime}\): \(\eta = q/p^{\prime}\). According to Wood et al. [72], the plastic deformation, which controls the softening response of soil, depends on the difference between current stress ratio and the peak stress ratio \(\eta_{\text{peak}}\) attainable at the current state. In order to determine the peak stress ratio and subsequently address the softening issue of dense sand in the drained condition, Wood et al. [72] proposed a relation that links the peak stress ratio to the state parameter: \(\eta_{\text{peak}} = M - n\psi_{\text{peak}}\), which was later modified by Manzari and Dafalias [41] as \(\eta_{\text{peak}} = Me^{{ - n\psi_{\text{peak}} }}\). This relation indicates that when the material is in a dense state (\(\psi < 0\)), \(\eta_{\text{peak}}\) is higher than \(M\) and its value depends on the parameter \(n\). In this study, we use the relation \(\eta_{\text{peak}} = \alpha Me^{{ - n\psi_{\text{peak}} }}\), which better fits our numerical data. Simulation results for medium dense samples, shown in Fig. 8, indicate that \(\alpha\) is close to 1 and this relation is not controlled by the fines content but by the particle rolling resistance or particle shape in other words. In contrast, the parameter \(n\) widely alters (i.e. \(n\) decreases from 5.01 in Case A to 3.48 in Case C) in response to the change in rolling resistance of coarse and fine particles. This trend suggests that when the rolling resistance of both coarse and fines particles in the mixture increases, the stress parameter at peak reduces and vice versa. The reason is that at low particle angularity, the mixture has a narrower range of void ratio as they reach the loosest state at lower void ratio [8]. If two samples have the same state parameter, e.g. the same value of \(\psi\), and are under the same loading path, the sample with lower rolling resistance has higher relative density compared to the other [22] due to less capacity to dilate. Therefore, the stress ratio \(\eta\) at peak of this sample with a lower rolling resistance is also higher, while its value of \(M\) at critical state is smaller due to less contribution from rolling resistance (see results in Sect. 3.1). Therefore, the stress parameter at peak \(\left( {\eta_{\text{peak}} /M} \right)\) dramatically increases when particles in the mixture have low rolling resistance or are more rounded.

Fig. 8
figure 8

Rolling resistance effect on the stress parameter at peak state of silty sand

Comparing results among Cases B, C, and D, as shown in Fig. 8, relatively increasing the rolling resistance of fines particles compared to that of coarse particles also reduces the peak stress parameter \(\left( {\eta_{\text{peak}} /M} \right)\). In these three cases, coarse particles are kept to the same anti-rotation state (see Table 1), and rolling resistance of fines in Case C is smaller than that in Case D and larger than in Case B. As a result, the value of \(n\) in Case C \(\left( {n = 3.49} \right)\) is smaller than that in Case B \(\left( {n = 3.69} \right)\) and higher than that in Case D \(\left( {n = 3.16} \right)\). This is due to the mutual effect between volume and shear responses on the overall behaviour at the peak state. Shape or anti-rotation effect alters the shear strength. At the same time, it alters the density state at the same state parameter because of the changes in the maximum and minimum void ratio. Therefore, in these three cases, the variation of parameter \(n\) is not as much as in Case A. Actually, if the percentage of fines increase up to 100%, \(n\) in Case B and Case A are equivalent. This implies that the relation between \(\eta_{\text{peak}} and \psi_{\text{peak}}\) is not completely independent of fines percentage. However, an approximation can be applied for silty sands with low fines content. A similar point of view can be found in the study of Wei and Yang [69].

The peak state is important to determine the hardening rate of the material because, as at this state, the hardening modulus is considered to be zero. Moreover, the hardening rate must be proportional to the state parameter to the peak not to the critical state [28]. Therefore, it is pivotal to determine the state parameter at peak, \(\psi_{\text{peak}}\). Figure 9 shows that the state parameter at peak can be calculated from the initial state parameter \(\psi_{0}\) via a linear relationship: \(\psi_{\text{peak}} = \gamma \psi_{o}\) when \(\psi_{0} < 0\) or samples are dense and \(\gamma < 1\). The values of \(\gamma\) are 0.87, 0.74, 0.70, and 0.68 for Cases A, B, C, and D in Fig. 9 correspondingly. The variation of \(\gamma\) is similar to the \(\left( {\eta /M} \right)_{\text{peak}}\) and can be explained similarly by the complex mutual effect between volume and shear responses on the overall behaviour at the peak state.

Fig. 9
figure 9

The linear relationship between the initial state parameter and state parameter at peak

Overall, the peak state strongly depends on the particle rolling resistance or the particle shape. The particle shape represented by the particle rolling resistance causes a similar impact on the change of \(\psi_{\text{peak}}\) and \(\left( {\eta /M} \right)_{\text{peak}}\). The simulation results suggest that with a certain initial condition and a known critical state, we can determine the locus of peak state and the shear strength at peak, which can be considered independent of fines at the low fines content. This could be applied to improve the constitutive model for a more accurate prediction of the behaviour of silty sands.

Dilatancy

It is well known that failures can occur progressively, with different mobilisations of strength and dilatancy at various locations along a developing slip surface. Therefore, dilatancy towards critical states is central to an understanding of soil behaviour. Both effective stress and soil density affect the rate of dilatancy of soils and thereby their strength parameters. And as previously discussed, the maximum rate of dilation is usually associated with the peak strength. Thus, the maximum dilatancy angle \(\theta_{\hbox{max} }\) can be calculated by the maximum ratio of axial-strain increment \({\text{d}}\varepsilon_{1}\) and volumetric strain increment \({\text{d}}\varepsilon_{v}\) [67]:

$$\sin \left( {\theta_{\hbox{max} } } \right) = \frac{{\begin{array}{*{20}c} {\left( {\frac{{{\text{d}}\varepsilon_{v} }}{{{\text{d}}\varepsilon_{1} }}} \right)_{\hbox{max} } } \\ \end{array} }}{{2 - \left( {\frac{{{\text{d}}\varepsilon_{v} }}{{{\text{d}}\varepsilon_{1} }}} \right)_{\hbox{max} } }}$$
(4)

The maximum dilatancy angle in relation to confining pressure and fines content is plotted in Fig. 10. The rolling resistance of fines particles causes a similar effect on the dilatancy of host sand as it does on the critical state. Adding fines with higher rolling resistance will increase the maximum dilatancy angle and vice versa. This is only applicable when fines content is smaller than a threshold fine content which is close to TFC, because, within this range, fines have more positive attribution on the structure of the mixture, which is mainly determined by the sand matrix. This agrees well with the experimental results reported by Yang and Wei [66], Thevanayagam and Shenthan [42] and Xiao et al. [76]. The results also show that the maximum dilatancy angle increases with the decrease in confining pressure regardless of the amount of fines or fines rolling resistance. However, the slope of the increment is strongly dependent on the amount of fines added and its relative rolling resistance to coarse. If rolling resistance of fines particles is higher than that of coarse particles, the maximum dilatancy angle rapidly increases when increasing fines percentage up to the threshold fines content. For a higher amount of fines added after this transition fine content, the dilatancy of the mixture reduces or remains unchanged depending on the initial confining pressure. If fines rolling resistance is higher than coarse rolling resistance, when increasing the amount of fines, the maximum dilatancy angle remains unchanged at low initial confining pressure. Otherwise, if fines rolling resistance is equivalent to or smaller than coarse rolling resistance, adding fine will lower the dilatancy effect. The threshold fine content, where the attribution of fines ceases the effect on the maximum dilatancy angle, is also changed with the rolling resistance difference between coarse and fine particles. Figure 10 shows that higher rolling resistance of fines will lower the value of transition fines content, besides the effects of confining pressure on dilatancy. About 5 degrees decrease in maximum dilatancy angle if confining pressure increases from 100 to 500 kPa. The increment is almost similar for different rolling resistances of fines.

Fig. 10
figure 10

Max dilatancy angle with different fines’ rolling resistances: a comparison among 3 different rolling resistances of fines for Cases B, C, and D; b Case D; c Case C; and d Case B

The experimental observations reported by Oda and Kazama [27] strongly suggest that rolling, rather than sliding, is the dominant micro-deformation mechanism controlling the dilatancy. Dilatancy parameter (\(D\)) is often considered as a function of both internal and external states of granular materials. When subjected to shear, loose sample contracts, and dense sample dilates, and according to critical state soil mechanics, a loose or dense state is not only defined in terms of density but also of the confining pressure. Therefore, Li and Dafalias [38] stated that dilatancy must be in general controlled by stress ratio \(\eta\), state parameter \(\psi,\) and intrinsic material properties. They proposed the following dilatancy relation:

$$D = D_{o} \left( {e^{{A_{d} \psi }} - \frac{\eta }{M}} \right)$$
(5)

where \(D_{o}\) and \(A_{d}\) are constants, which control the dilatancy behaviour of soil. Experimentally, the \(A_{d}\) can be inferred by data fitting method at the zero-dilatancy state:

$$e^{{A_{d} \psi_{D = 0} }} = \frac{{\eta_{D = 0 } }}{M}$$
(6)

Therefore, zero-dilatancy state is one of the most important states reflecting the shearing response of soils. At zero-dilatancy state, also characterised as the characteristic state, there is temporarily no volumetric deformation. This is a transition from elastoplastic contraction at the beginning of shear to the dilation behaviour happening along the path to the critical state on the initial medium-to-dense soil. Figure 11 shows that there is a linear relationship between the initial state parameter and the characteristic state parameter as follows:

$$\psi_{D = 0} = A\psi_{0} + B$$
(7)
Fig. 11
figure 11

The linear relationship between the initial state parameter and state parameter at zero dilatancy

The level of rolling resistance of fines and coarse contents is less likely to affect the relative position of the transition compared to the initial value of the state parameter. Comparing different rolling resistances of fines and coarse contents, and different amounts of fines added as shown in Fig. 11, it can be seen that values of \(A\) and \(B\) in Eq. (7) somewhat decrease when increasing the rolling resistance of the mixture. The negative value of B implies that if the initial state is the critical state or state parameter is zero, the sample still contracts before going back to the critical state.

The DEM simulations suggest that dilatancy also depends on the stress state. The stress state is represented in Eq. (5) via stress parameter \(\eta /M\). To consider the effect of rolling resistance on the stress parameter at zero-dilatancy state, rather than Eq. (6), we use a more general equation \(\left( {\eta_{D = 0} /M} \right) = \beta e^{{A_{d} \psi_{D = 0} }}\). The exponential relationships between stress ratio and state parameter for different fines contents and fine rolling resistances at zero dilatancy are plotted in Fig. 12. These plots do not vigorously support Eq. (6) especially in a dense state. The denser the sample is, the more divergent the stress parameter becomes. This implies that Eq. (6) cannot unify the dilatancy of the granular mixture with varying fines content. A similar observation has also been reported in the experimental study of Wei and Yang [69]. They ascribed this non-unification to the particle shape and gradation. This observation is agreed well with our simulation results that increasing rolling resistance of fine and coarse particle will reduce the data scatter in the correlation between stress ratio and state parameter at zero-dilatancy state. Moreover, the simulation results can likeably demonstrate the effect of rolling resistance on the dilatancy in the way that enhancing the rolling resistance of both coarse and fines will reduce the stress ratio at zero dilatancy.

Fig. 12
figure 12

Rolling resistance effect on the zero-dilatancy state of silty sand

The above results suggest that the particle angularity (or rolling resistance) strongly affects the stress ratio at zero-dilatancy state. Samples containing rounded particles tend to have lower elastoplastic contraction potential. Otherwise, samples containing more angular particles (or high rolling resistance) have more enduring structure due to the stronger stick–slip contacts among particles. That also explains why samples with high angular particle reach the zero-dilatancy state at a higher volumetric deformation. However, the scatter data with value \(\beta\) far from 1 shown in Fig. 12 imply that the fines content, or the particle size distribution in general may overweight the particle shape effect in controlling the zero-dilatancy state. This requires a better relationship to unify dilatancy behaviour with various fines contents. Nevertheless, the fines content is less likely to control the state parameter at zero dilatancy, which is important in predicting the behaviour of silty sand especially in cyclic loading problems where the initial state changes after each loading cycle.

Effect of particle rolling resistance on undrained behaviour of silty sands

Typical undrained behaviour of sand can be found in Fig. 13 (Case A). The undrained behaviour can be classified into three different types: non-flow, limited flow, and flow. Similar to drained behaviour, undrained behaviour highly depends on the initial void ratio of sand. Thus, non-flow, limited flow, and flow responses are associated with dense sand, medium dense sand, and loose sand, respectively. Loose sand exhibits the flow behaviour where effective shear stress quickly reaches the peak value at small shear strain and then softens until reaching the equilibrium at a critical state, also known as the steady state in undrained tests. For medium dense sand, after initial softening, hardening responses take place until reaching its steady state. There is no softening occurring in shearing of dense samples; instead, dense samples continuously dilate until reaching their steady state. Figure 13 shows that fines particles of silt size can significantly alter the flow potential of silty sands. In general, at similar initial void ratios, adding fines up to TFC will continuously increase the liquefaction potential of the sands. However, further increasing fines content after TFC will reduce the possibility of liquefaction. In addition, the effect of fines on the liquefaction potential depends on not only the added amount but also the anti-rotation magnitude of fines particles. This observation agrees well with the experimental results reported in the studies of Thevanayagam et al. [62], Polito and Martin [52], Xenaki and Athanasopoulos [74], Yang et al. [81], and Papadopoulou and Tika [48]. Comparing results for Case A and Case C shown in Fig. 13, it can be concluded that the sand–silt mixture with well-rounded particles has the higher collapsibility in response to the undrained shear than the one with higher angular shape for both coarse and fines particles.

Fig. 13
figure 13

Evolution of undrained shear stresses and pore water pressure for different simulation cases (void ratio e ≈ 0.55 for Case A, and e ≈ 0.58 for Cases B, C , and D)

On the other hand, comparing among Cases B, C, and D in Fig. 13 where coarse particles are set at the same rolling resistance shows that shear strength increases when increasing fines with higher rolling resistance. At the same percentage of fines and the same initial void ratio, the silty sand which has higher anti-rotation fines exhibits a reduced-strain softening response and a decrease in liquefaction potential. This strong correlation between particle rolling resistance and mechanical behaviour explains the very diverse views as to whether the effect of fines is negative or positive on the liquefaction resistance of silty sands ([36, 45, 61]). The detailed effect of particle rolling resistance on the undrained steady state, undrained peak state, and phase transformation state will be discussed in the next sections.

Undrained critical state

There are contradictory conclusions found in the literature about the correlation between critical state and drainage condition. On the one hand, it is assumed that the critical state for soil that has the same intrinsic properties is uniquely irrespective of drainage conditions. This comment associates closely with many experimental studies conducted by Sladen et al. [58], Been et al. [2], Bobei et al. [3]. On the other hand, Yamamuro and Lade [36] found that the slope of the drained critical state was higher than that of the undrained critical state in triaxial compression tests, especially at low confining pressure.

Simulation results in this study shown in Figs. 14 and 15 support the viewpoint that the critical state value in the stress space is almost the same for both drained and undrained loading conditions. It is perspicuously shown in Fig. 14a that a large difference between the anti-rolling of coarse and fines particle leads to a considerable change of the critical state slope \(M\) with altering fines content. However, there is a small difference in the critical state slope amplitude between two different drainage conditions. The maximum critical void ratio \(e_{\varGamma }\) also shows a similar response to drained and undrained loading conditions when its value hardly varies with drainage, as shown in Fig. 14b.

Fig. 14
figure 14

Drainage condition effect on the coarse–fines mixture with different rolling resistances of coarse and fines particle: a critical state slope in the shear plane; b maximum critical void ratio

Fig. 15
figure 15

Effect of drainage condition in the slope of critical state in volumetric plane

Nevertheless, the slope of the critical state line \(\lambda_{c}\) in the volumetric plane is slightly affected by the drainage condition and the rolling resistance of fines added into the mixture. At the TFC around 20% of fines, the drained \(\lambda_{c}\) is closest to the undrained counterpart. TFC forms the optimal particle distribution where fines particles are efficiently filling the pore space between coarse particles; therefore, shearing with or without volume constraints is less different from each other. However, if the fines content is far from the TFC, increasing the rolling resistance of fines reduces the rotation of the volumetric critical state line in undrained condition evidenced in the comparison among Case B, Case C, and Case D in Fig. 15. Although the reduction is not very significant, it is showing that the volume constraint enhances the influence of particle shape on the volume change resistance of the sands. This is because the constraint in volumetric deformation regulates the motion of particles. Although there is space when pores are not fulfilled by fines when fine content \(f\) is smaller than TFC, particles cannot freely occupy that space because their motion is restricted following the undrained condition. As a result, higher confining pressure is required to reach the same void ratio in undrained condition compare to drained one. Inversely, if FC is higher than TFC, fines motions govern the material flow as fines dominate the soil structure. The less constraint in the motion leads to a stronger change in the soil structure. Therefore, the flow of fines among coarse particle is stronger in drained condition compared to the undrained. At this time, coarse particles act as the obstacles hindering the fines flow. As a result, the compressibility of soil in the drained condition is lower than in undrained condition.

Phase transformation state

In undrained tests, medium loose to dense sand often exhibits contraction followed by dilation during shearing, which is termed as phase transformation behaviour. A phase transformation line can be defined as a line passing through phase transformation states of various undrained conditions as in Fig. 16a. Different from drained tests, dilatancy at the transition point in undrained condition cannot be characterised by \(d\varepsilon_{v} = 0\) as this is a constraint required for any undrained tests. To define the phase transformation point in undrained tests, the change in pore water pressure (\(\Delta u\)), which is analogous to the change in \(\varepsilon_{v}\) in drained test [46, 60], is used. At this point, the change in pore water pressure vanishes, \(\Delta u = 0\).

Fig. 16
figure 16

Effect of rolling resistance on phase transformation lines: a phase transformation lines in Case D; b comparison among cases

The transition of phase transformation line under the effect of rolling resistance has the same tendency as that of characteristic line. In this respect, the phase transformation line rotates in the same direction as the critical state when changing the amount of the fine in the coarse–fines mixture. As shown in Fig. 16, higher rolling resistance of fines will rotate the phase transformation line counter-clockwise in \(p - q^{\prime}\) plane and vice versa. Otherwise, if fines particle has a similar degree of anti-rotation as coarse particles, adding more fines-to-silty sands hardly affect the locus of phase transformation state.

The phase transformation line plays a similar role to the undrained response as the characteristic line does for drained behaviour. The stress path will turn its direction after the phase transformation occurs. As shown in Fig. 13, at this point, the mean effective stress \(p^{\prime}\) reaches the minimum value. The increment of \(p^{\prime}\) turns to zero. Consequently, the elastic and plastic volumetric strain increments are also zero. Comparing the results shown in Figs. 13 and 6b, this study supports the idea that phase transformation state is comparable to the critical state [35].

Instability stress ratio and instability line

Instability state or pre-failure flow state refers to the rapid generation of large plastic strain or large excess pore water pressure and a dramatically decreased shear strength under undrained strain-controlled triaxial tests [82]. This state has been considered as one of the failure mechanisms that lead to the collapse of granular soil slopes in a number of case studies [82]. For specimens exhibiting strain-softening behaviour in undrained condition, an instability line can be defined by connecting the peak stress obtained from a series of effective stress paths in undrained tests. As shown in Fig. 17, the instability line is not unique. It varies with the void ratio of sand and is affected by the effective stresses in general [9, 82]. Therefore, it is challenging to find a clear pattern in the relationship between the initial state parameter and the stress ratio at the instability state when altering the fines content. However, it can be seen from each case in Fig. 17 that adding fines up to a certain value tends to increase \(\eta_{\text{instability}}\); further increasing fines after this value, \(\eta_{\text{instability}}\) will decrease. Moreover, this value increases with the increase in the rolling resistance of the mixture. Comparing among cases in Fig. 17, the effect of the rolling resistance on the instability state of the mixture can be seen in the transition of the \(\psi_{0} - \eta_{\text{instability}}\) relation. For the same fines content, increasing rolling resistance of coarse particles, fines particles, or both will shift \({\text{the}}\; \psi_{0} - \eta_{\text{instability}}\) line to the right. The instability state is reached at looser state or higher stress ratio. This implies that the mixture is more stable as rolling resistance increases.

Fig. 17
figure 17

Stress ratio at instability state—straight lines represent the interpolation curves \(\eta_{\text{instability}}\)

The instability line may be used to identify the stress states in the \(p^{\prime} - q\) stress plane beyond which the instability behaviour of the granular soils commences. The region between the instability line and the critical state line has been called the zone of potential instability [34]. The behaviour of granular material is highly sensitive with the void ratio in undrained triaxial simulation and is more sensitive in the potential instability zone. Because most of the constitutive models for soils do not often focus on dilatancy behaviour at the instability state, predicting the peak value at instability is not precise. Locus of instability line is normally over-predicted as it happens sooner than observed in experimental and numerical simulations [54, 84]. This, therefore, become a big accumulative error in cyclic loading. Our simulation results with different rolling resistances shown in Fig. 18 point out that although there is no unique relationship between stress ratio and state parameter at the instability point, there is a linear correlation between initial state parameter and state parameter at instability expressed as \(\psi_{\text{instability}} = \psi_{0} + l\). In another words, instability state is potentially reached at a constant distance from initial state and changing fines content or particle size distribution hardly influences this volumetric potential. This constant distance \(l\) is, however, affected by the rolling and sliding resistances of particle, which can either strengthen or weaken the inter-particle contact. Increasing the rolling resistance of mixture will increase the value of \(l\) and vice versa.

Fig. 18
figure 18

State parameter at instability state

As the instability state refers to elastic-to-elastoplastic transition, which takes place before phase transformation state and before the soil structure significantly change under shear, it is important to predict well this state of the material. Our simulation results suggest that for a better prediction of liquefaction behaviour of soils, constitutive models should include the distance to the instability state and update this distance every cycle of cyclic loading.

Conclusions

The effect of rolling resistance on the mechanical behaviour of silty sand was investigated to provide a better understanding of how particle properties relating to particle shape alter the macro-behaviour of the material. A series of triaxial simulations with different rolling resistances for fines and coarse particles was conducted for sand–silt mixture with varying fines contents in different initial states and drainage conditions. The salient conclusions that can be drawn from this work are as follows:

  1. 1.

    The elastic–plastic spring–dashpot rolling resistance model can be employed to adequately approximate the influence of particle shape on the mechanical behaviour of granular materials. The DEM simulations with the rolling resistance model successfully capture the state-dependent behaviour of granular soil, including contraction and dilation, instability and quasi-steady state in medium dense soil, and dilatancy and hardening in dense soil with different initial densities, confining pressures and loading conditions.

  2. 2.

    The rolling resistance (or particle shape) has a small effect on the transition fines content, the percentage of coarse–coarse, coarse–fines, and fines–fines contacts in the mixtures, and the rotation of the critical state in the volumetric plane. The critical state lines with different fines contents are almost parallel to each other in the volumetric plane regardless of rolling resistance. However, the rolling resistance strongly affects the maximum void ratio at the critical state. The critical state in the stress plane is independent of drainage condition. Nevertheless, the results show that the rotation of critical state in the volumetric plane is less significant in undrained than in drained conditions when fines content varies.

  3. 3.

    Rolling resistance can enhance or reduce the shear strength of sand–silt mixture depending on the ratio between the rolling resistances of fines and coarse particles. Fines with rolling resistance higher than that of coarse contents will enhance the shearing resistance of the host sand as it increases both the critical state slope in stress space and peak friction angle. In contrast, fines with rolling resistance with lower than or similar to that of coarse content will reduce the shear strength of the mixture or have an insignificant effect on the mechanical behaviour of the host sand.

  4. 4.

    The correlation between peak stress ratio and state parameter at peak is almost independent of fines content. However, this correlation strongly depends on the rolling resistance of fines and exhibits a similar trend to the critical state when changing the rolling resistance of fines. The characteristic state in the drained triaxial test and the phase transformation state in the undrained triaxial test, which are similar to each other, also show a similar response to the fines rolling resistance.

  5. 5.

    Rolling resistance has a similar effect on the phase transformation state as it does on the critical state in the stress plane. Although there is no clear effect of rolling resistance on the instability line, it has been shown that there is a constant distance between the initial state parameters and the instability state. This correlation is independent of the anti-rotation effect and can be further applied to improve constitutive models in predicting the monotonic and cyclic undrained behaviour of soil.

  6. 6.

    Liquefaction resistance potential increases with increasing fines content up to transition fines content and decreases thereafter with further increasing fines content. This trend is independent of the shape of particles of the mixture as the undrained behaviour strongly coincides with the volume response of soils. However, comparing with adding fines of a more angular shape, adding the same amount of more rounded fines leads to higher collapsibility or liquefaction susceptibility. The results shown in this study confirm that grain shape represented by rolling resistance is the cause of experimentally observed different shear behaviours in static liquefaction of silty sands.

Finally, the paper also demonstrates that DEM with rolling resistance could provide consistent results and a general understanding of the effect of particle shape on the behaviour of silty sand. However, the employment of rolling resistance for spherical particles still has limitation compared to using real particle shapes. One approach to mitigate this adverse impact is to replace coarse particles by particle clumps, while fine particles can still be modelled using rolling resistance. This approach will be further explored in future works.